Simulation study of earthquakes based on the two-dimensional Burridge-Knopoff 

model with the long-range interaction 
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Spatiotemporal correlations of the two-dimensional spring-block (Burridge-Knopoff) models of 
earthquakes with the long-range inter-block interactions are extensively studied by means of numer- 
ical computer simulations. The long-range interaction derived from an elastic theory, which takes 
account of the effect of the elastic body adjacent to the fault plane, falls off with distance r as 1/r 3 . 
Comparison is made with the properties of the corresponding short-range models studied earlier. 
Seismic spatiotemporal correlations of the long-range models generally tend to be weaker than those 
of the short-range models. The magnitude distribution exhibits a "near-critical" behavior, i.e., a 
power-law-like behavior close to the Gutenberg-Richter law, for a wide parameter range with its 
B-value, B ~ 0.55, insensitive to the model parameters, in sharp contrast to that of the 2D short- 
range model and those of the ID short-range and long-range models where such a "near-critical" 
behavior is realized only by fine-tuning the model parameters. In contrast to the short-range case, 
the mean stress-drop at a seismic event of the long-range model is nearly independent of its magni- 
tude, consistently with the observation. Large events often accompany foreshocks together with a 
doughnut-like quiescence as their precursors, while they hardly accompany aftershocks with almost 
negligible seismic correlations observed after the mainshock. 

PACS numbers: 91.30.Px,05.10.-a 
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I. INTRODUCTION 



An earthquake is a stick-slip dynamical instability of a 
pre-existing fault driven by the motion of a tectonic plate 
[H> 0- While an earthquake is a complex phenomenon, 
certain empirical laws such as the Gutenberg-Richter 
(GR) law and the Omori law concerning its statistical 
properties are known to hold. Understanding the origin 
of such statistical properties of earthquakes is one of im- 
portant issues left in earthquake studies. As a useful tool 
in such studies, many researchers have used the so-called 
spring-block model originally proposed by Burridge and 
Knopoff (BK) In this model, an earthquake fault is 
simulated by an assembly of blocks, each of which is con- 
nected via the elastic springs to the neighboring blocks 
and to the moving plate. All blocks are subject to the 
friction force, the source of the nonlinearity in the model, 
which eventually realizes an earthquake-like frictional in- 
stability. While the spring-block model is obviously a 
crude model to represent a real earthquake fault, its sim- 
plicity enables one to study its statistical properties with 
high precision. 

Carlson, Langer and others 0, [E S S B studied 
the statistical properties of the ID and 2D BK mod- 
els quite extensively, paying particular attention to the 
magnitude distribution of earthquake events. The spring- 
block model has also been extended in several ways, e.g., 
taking account of the effect of viscosi ty JlO. [ill. Il2l| , mod- 
ifying the form of the friction force [id, Jl2j. Il3j|. driving 
the system only at one end of the system [14j , or by incor- 
porating the rate- and state-dependent friction law [TH ]. 
The present authors studied in the previous papers the 
statistical properties of the ID and 2D BK models, fo- 
cusing on their spatiotemporal correlations (l6l . [TtI [l8| . 



These studies have revealed several interesting features 
of the ID and 2D BK models. 

Meanwhile, the BK models studied in most of the 
previous works assumed that the inter-block interaction 
works only between nearest-neighboring blocks. This cor- 
responds to the situation where a thin isolated plate is 
subject to the friction force and is driven by shear force 
jl9j. However, a real fault is not necessarily a thin iso- 
lated plate, and the elastic body extends in a direction 
away from the fault plane. Considering the effect of 
such an extended elastic body adjacent to the fault plane 
amounts to considering the effective inter-block interac- 
tion to be long-ranged. In order to make the model more 
realistic, it is important to take account of effect of the 
long-range interaction, together with the effect of the di- 
mensionality of the fault. In this connection, we note 
that, in the study of thermodynamic phase transition in 
equilibrium, it has been wellknown that the spatial di- 
mensionality and the range of the interaction are major 
elements affecting the universality class of the transition. 

Hence, in the present paper, we study the statisti- 
cal properties of the 2D BK model with the long-range 
inter-block interaction derived from an clastic theory, in 
comparison with those of the BK models with the short- 
range (nearest-neighbor) interactions studied earlier, in 
order to get information how the long-range nature of 
the interaction, expected to arise from the elastic prop- 
erties of the crust adjacent to the fault plane, affects the 
statistical properties of earthquakes. 

We assume that the 3D elastic body, where the 2D BK 
models with the long-range interaction is supposed to lie, 
are isotropic, homogeneous and infinite. A fault surface 
is assumed to be a plane lying in this elastic body and to 
slip along one direction only. As a further simplification, 
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we adopt a static approximation for an elastic equation 
of motion describing the elastic body. This assumption 
is justified when the velocity of the seismic-wave propa- 
gation is high enough compared with the velocity of the 
seismic-rupture propagation. As shown in the appendix, 
these assumptions give rise to a spring constant between 
blocks decaying with their distance r as 1/r 3 . 

Certain properties of the BK model with the long- 
range interaction, or the BK model extended in the direc- 
tion orthogonal to the fault plane, were already studied. 
These include the 2D BK model extended in the direc- 
tion orthogonal to the fault plane [2(J, the 2D cellular 
automaton version of the BK model with the long-range 
interaction decaying as 1 /r 3 [2l[ . In particular, Xia et al 
recently studied the ID BK model with a variable range 
interaction where a block is connected to its R neigh- 
bors with a rescaled spring constant proportional to 1/R 
[22L I23I ] . The type of the long-range model considered by 
Xia et al may be regarded as a mean-field type, since the 
model reduces to the mean-field infinite-range model in 
the R — > 00 limit. 

In the present paper, we extend our previous studies 
on the spatiotemporal correlation properties of the short- 
range BK models [HI, [l?], HI], we investigate the spa- 
tiotemporal correlation properties of the 2D BK model 
with the long-range power-law interaction derived from 
an elastic theory which is expected to capture the effect 
of the elastic body adjacent to the fault plane. Our work 
can also be regarded as an extension of the recent work of 
Xia et al [22I |23| : First, we extend the model dimension- 
ality from ID to more realistic 2D. Second, we consider 
the long-range interaction derived from an elastic theory, 
decaying as a power law with distance, which is different 
from the mean-field-type long-range interaction consid- 
ered in Ref . [22I [23[. Third, we calculate various Spa- 
tiotemporal correlation functions to further examine the 
properties of seismicity under the influence of the long- 
range interaction. In view of the situation that many of 
the previous works on the BK model were performed for 
the ID model, however, we also perform for comparison 
a similar numerical analysis complementally for the ID 
BK model with the long-range power-law interaction. 

The present paper is organized as follows. In §11, we 
introduce the model and explain some of the details of 
our numerical simulation. The results of our simulations 
on the 2D BK model with the long-range interactions are 
presented in §111. We show the results of the event-size 
distribution, the mean displacement, the mean number 
of failed blocks and the mean stress drop at a seismic 
event, together with various types of spatiotemporal cor- 
relation functions of seismic events, including the local 
recurrence-time distribution, the seismic time-correlation 
function before and after the mainshock, the time devel- 
opment of the seismic space-correlation function before 
and after the mainshock, and the time development of the 
magnitude distribution function before the mainshock. 
The derivation of the long-range inter-block interaction 
from an elastic theory is given in Appendix A. The re- 



sults of our calculation on the ID BK model with the 
long-range power-law interaction is also presented in the 
Appendix B. Finally, §IV is devoted to summary and 
discussion. 



II. THE MODEL AND THE METHOD 

First, we describe the 2D BK model with the nearest- 
neighbor interaction. The 2D BK model represents a 
"fault plane" by an assembly of blocks, which is taken 
to be an x — z plane consisting of a 2D square array 
of blocks containing N x blocks in the x-direction and 
N z blocks in the z-direction. All Blocks are assumed to 
move only in the ^-direction along strike, and are subject 
to the friction force $. Each block is connected with its 
four nearest-neighbor blocks via the springs of the elastic 
constant k c , and is also connected to the moving plate 
via the spring of the elastic constant k p . 

In the simplest case where the interaction works 
only between the nearest-neighbor blocks in a spatially 
isotropic manner, the equation of motion of the block at 
site is given by 

mUij = k v (v't' - Ui.j) + k c (U i+ i t j + Ui ij+ i 

+Ui-u + Uij-! - 4U id ) - §(U id ), [ 1 

where m is the mass of a block, t' is the time, Uij is the 
displacement along the ^-direction of the block at site 
and v' is the loading rate representing the speed 
of the plate. The equation is made dimensionless in the 
same way as in [17], i.e., the time t' is measured in units 
of the characteristic frequency u> — ^Jk p /m and the dis- 
placement Ui.j in units of the length L = $>(0)/k p , $(0) 
being a static friction. Then, the equation of motion can 
be written in the dimensionless form as 

iii = vt- Uij + l 2 (u l+1J + u hj+1 
+u i -ij+Ui i j- 1 -4ui 1 j)-</)(ui), 

where t — t'ui is the dimensionless time, Uij = Uij/L 
is the dimensionless displacement of the block 
I = \Jk c /k p is the dimensionless stiffness parameter, 
v = v'/(Lu>) is the dimensionless loading rate, and 
4>(ui) = 3>(f7i)/$(0) is the dimensionless friction force. 

The nearest-neighbor model mentioned above neglects 
the effect of the elastic body in a direction away from the 
fault. As shown in Appendix A, taking account of this 
effect amounts to taking the inter-block interaction to be 
long-ranged. The interaction between the two blocks at 
sites (i,j) and is given in the dimensionless form 

by 

which falls off with distance rasl/r 3 . Then, the equation 
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of motion of the 2D long-range can be written as 



= v t — Ui 



2 \5'-0\' 



i) 



(4) 

If one restricts the range of interaction to nearest neigh- 
bors and takes the spatially anisotropic spring constant 
to be isotropic, l x = l z = I, one recovers the isotropic 
nearest- neighbor model described by Eq.(2). 

The "isotropy" assumption l x = l z is equivalent to 
putting the Lame's constant to vanish, A = 0. In fact, 
the possible effect of such spatial anisotropy of the 2D 
BK model was studied within the nearest-neighbor in- 
teraction in our previous paper [lsj . It was observed 
that the property of the anisotropic model was close to 
the corresponding isotropic model characterized by the 
mean spring constant I = (l x + l z )/2 so that the spa- 
tial anisotropy did not cause any qualitative new fea- 
ture on the statistical properties of the model. Thus, 
in the present paper, we put l x = l z = I for simplicity. 
The investigation of the recurrence-time distribution of 
the anisotropic model with l x ^ l z was recently made in 
Ref. within the nearest-neighbor interaction. 

In the present paper, we also discuss in the appendix 
the properties of the ID BK model with the long-range 
interaction, to clarify the role of the model dimensional- 
ity and to make comparison with the previous works on 
the various ID BK models. We derive the ID BK model 
with the long-range interaction from the corresponding 
2D model by imposing the constraint that the systems is 
completely rigid along the ^-direction corresponding to 
the depth direction, i.e., u(x, z, t) — u(x, t). As shown in 
Appendix A, this yields an effective inter-block interac- 
tion decaying with distance r as 1/r 2 , 
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Then, the equation of motion of the ID BK model may 
be given in the dimensionless form by 



Hi = vt 



i'^i \i — V | 2 



(6) 



As the form of the friction force 4>, we use a simple 
velocity-weakening friction force which is a single- valued 
function of the velocity. As its explicit functional form, 
we use the form introduced by Carlson and Langer Q, 



(~oo,l], 

l-CT 

l+2ail/(l-u) • 



for ii < 0, 
for ii > 0, 



(7) 



where the friction force immediately drops to 1 — a on 
sliding, and decays toward zero with a rate proportional 
to the parameter a as the velocity increases. The back- 
slip is inhibited by imposing an infinitely large friction for 
iii < 0, i.e., <j)(u < 0) = — oo. This friction force repre- 
sents the velocity-weakening friction force. Although real 
friction force is of course more complex, not depending 



on the velocity alone [l[ , we use the friction force (7) for 
simplicity. 

The friction force is characterized by the two parame- 
(itij^rs, a and a. The former, a, represents an instantaneous 
drop of the friction force at the onset of the slip, while the 
latter, a, represents the rate of the friction force getting 
weaker on increasing the sliding velocity. The a = case 
represents the simplest Coulomb friction law where the 
friction force instantaneously drops from the static value 
1 to its dynamical value 1 — a as soon as the block begins 
to slide, and is kept constant on sliding irrespective of 
the velocity. The a = oo case also corresponds to the 
another Coulomb friction law where the dynamical fric- 
tion immediately drops to zero on sliding. In addition 
to these frictional parameters, the model possesses one 
more material parameter, an elastic parameter I. 

In the present paper, we try to cover a rather wide 
range of the parameter a in the range a = [0, oo], and 
systematically examine the a-dependence of the results. 

We also assume the loading rate v to be infinitesimally 
small, and put v = during an earthquake event, a very 
good approximation for real faults [6|. Taking this limit 
ensures that the interval time during successive earth- 
quake events can be measured in units of v~ x irrespective 
of particular values of v. 

A seismic event begins when the accumulated stress 
exceeds a static friction at one of the blocks in the sys- 
tem. Due to the effect of nonzero er, the block begins 
to move with a finite acceleration, which may (or may 
not) propagate to the neighboring blocks. The succes- 
sion of such propagating motion of blocks is regarded as 
a seismic event. The event is terminated when all blocks 
in the system come to rest again. The displacement of 
each block at an event is measured by the displacement 
of that block during the beginning and the end of this 
event. The condition of an infinitesimal a guarantees 
that no other event is triggered elsewhere in the system 
during the ongoing event. 

Numerical details are the same as in [ljj ■ We solve the 
equation of motion (4) or (6) by using the Runge-Kutta 
method of the fourth order, the width of the time dis- 
cretization At being At = 10~ 3 in most cases. The long- 
range interaction is summed over all blocks contained in 
the system. Total number of 10 5 ~ 10 7 events are gen- 
erated in each run, which are used to perform various 
averagings. The initial position of each block Ui(0) is 
generated randomly according to the uniform distribu- 
tion in the interval [0,0.02], with the zero initial veloc- 
ity Ui(0) ~ 0. In calculating the observables, initial 10 5 
events are discarded as transients. We judge whether 
the system reaches a stationary state by monitoring the 
stability of the magnitude distribution function (to be 
defined in detail below). 

In the 2D BK model, we follow [§[ and impose periodic 
boundary condition in the x-direction and free boundary 
condition in the z-direction, regarding the z-direction as 
the depth direction. For the most part of our calculation, 
the system size is taken to be N x = 160 and N z = 80 (or 
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N x = 60 and N z = 60). In the ID BK model studied in 
Appendix B, we impose periodic boundary condition. 



III. THE SIMULATION RESULTS 

In this section, we show the results of our numerical 
simulations on the 2D BK model with the long-range 
interaction for various observables. 



A. THE MAGNITUDE DISTRIBUTION 

We define the magnitude of an event of the 2D BK 
model, fi, as a logarithm of its moment M, 



Figs. 1(a) and (b) exhibit the computed magnitude dis- 
tribution function i?(/i) for smaller and larger values of 
a, i.e., (a) < a < 10 and (b) 10 < a < oo. The mag- 
nitude distribution R(fi)dfi represents the rate of events 
with their magnitudes in the range \pi, fi + dfj]. In the 
range a $ 1, only small events of /i IS 2 occur. As can 
be seen from Fig. 1(a), i?(/i) for smaller a 5 1 exhibits 
a near straight-line "near-critical" behavior over a cer- 
tain magnitude range, and drops off sharply at larger 
magnitudes. The associated B- value is estimated in the 
range 05a^ltobeB~ 0.59 from the slope of this 
straight line, which is rather insensitive to the change of 
the a-value. Of course, the observed behavior cannot be 
regarded as truly critical, since R(fJ-) drops off sharply 
beyond the threshold magnitude. 



= In M = In [ Aitij 



(8) 



where Ait^j is the total displacement during an event of 
the block at site (i, j) and the sum is taken over all blocks 
involved in the event. 
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FIG. 1: The magnitude distribution R(fi) of earthquake 
events of the 2D long-range BK model for the parameters 
I = 3 and a = 0.01. Fig. (a) represents R(fJ,) for smaller val- 
ues of the frictional parameter < a < 10, while Fig.(b) 
represents R(jJ,) for larger values of the frictional parameter 
10 < a < oo. The system size is 60 x 60. 
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FIG 2: The system-size dependence of the magnitude dis- 
tribution R(n) of earthquake events of the 2D long-range BK 
model for the parameters a — 30, I = 3 and a = 0.01. A 
sharp fall-off observed at larger magnitudes fi ?J 5 is not a 
finite-size effect. 

At a C 2, large earthquakes of their magnitudes fi ~ 8 
suddenly appear, while earthquakes of intermediate mag- 
nitudes, say, 2 5 /i 5 6, remain rather scarce. It means 
that large and small earthquakes are well separated at 
a ~ 2. Such a sudden appearance of large earthquakes 
at a = a c i ~ 2 coexisting with smaller ones has a feature 
of "discontinuous transition" . This feature is common 
to the case of the corresponding 2D short-range model 
[l8l |. On increasing a further, earthquakes of interme- 
diate magnitudes gradually increase their frequency. In 
the range of 2 a ^ 20, R(n) exhibits a "supercritical" 
behavior, i.e., exhibits a pronounced peak structure at 
a larger magnitude deviating from the GR law, while it 
still exhibits a near straight-line behavior corresponding 
to the GR law at smaller magnitudes. The existence of 
a distinct peak structure at a larger magnitude suggests 
that large earthquakes are more or less characteristic. 
Such a behavior of R(j-i) is sometimes called "supercriti- 
cal" , since R(fJ,) bends up at larger magnitudes (though 
it eventually falls off at still larger magnitudes). 
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FIG. 3: Comparison of the magnitude distribution R(u) of 
earthquake events of the 2D BK models with the long-range 
and the short-range interactions. Fig. (a) represents the case 
oi a — 10, while Fig.(b) the case of a = 30, with I = 3 and 
a — 0.01 being fixed. 



As a increases further, a characteristic peak becomes 
less pronounced and eventually vanishes at around a ~ 
25. R(p) exhibits again a near straight-line near-critical 
behavior over a wide magnitude range: See Fig. 1(b). At 
oe = a C 2 — 25, the associated £?-value estimated from 
the slope of this straight line is B ~ 0.55. The change 
from the supercritical to the near-critical behaviors at 
a = a c2 — 25 is continuous, in contrast to the discon- 
tinuous one observed at a — a c \ ~ 2. A very inter- 
esting observation here is that such a straight-line near- 
critical behavior persists even if a is further increased 
up to a = oo, and that the associated B-value is robust 
against the change of a. It should be noticed that this 
straight-line behavior of R(fi) cannot be regarded as a 
truly critical one, since R(fJ-) drops off sharply at very 
large magnitudes. This sharp fall-off of R(fi) observed at 
larger magnitudes fi ^ 5 is not a finite-size effect, as can 
clearly be seen from Fig. 2. 

Such a near-critical behavior realized over a wide pa- 
rameter range a c 25 is in sharp contrast to the behav- 
ior of the corresponding short-range model where 
at larger a > a C 2 exhibits a down-bending "subcritical" 
behavior, while a straight-line near-critical behavior is re- 
alized only by fine-tuning the a-value to a special value 



a ~ a C 2 ■ The robustness of the near-critical behavior of 
R(fji) observed in the 2D long-range model might have an 
important relevance to real seismicity. 




a 

FIG. 4: The phase diagram of the 2D BK models with the 
long-range and the short-range interactions in the frictional- 
parameter a versus elastic-parameter / plane. The parameter 
a is set to a = 0.01. To draw a phase diagram, the parameter 
range < a < oo and 1 < I < 10 is studied by simulations. 

In order to further illustrate the difference between the 
long-range and the short-range models, we compare in 
Fig. 3 R(fi) of the long-range and the short-range models. 
Fig. 3(a) represents the case of a = 10 in the supercrit- 
ical regime a c i < a < a C 2, while Fig. 3(b) represents 
the case of a = 30 in the near-critical regime a > a C 2- 
As can be seen from the figures, R(fj) of the long-range 
model exhibits much more pronounced straight-line be- 
havior mimicking the GR-law over a wider magnitude 
range, as compared with R(fJ-) of the short-range model. 

In Fig. 4, we summarize the behavior of in the 

form of a "phase diagram" in the frictional-parameter 
a versus the elastic-parameter I plane for the case of 
a = 0.01. As can be seen from the figure, the phase di- 
agram of the long-range model consists of three distinct 
regimes, two of which are near-critical regimes and one is 
a supercritical regime. The "phase boundary" between 
the smaller-a near-critical regime and the supercritical 
regime represents a "discontinuous transition" , while the 
one between the larger-a near-critical regime and the su- 
percritical regime represents a "continuous transition" . 
The "transition" between these different "phases", i.e., a 
near-critical phase for small a, a supercritical phase for 
intermediate a, and another near-critical phase for large 
a, is primarily dictated by the a-value. Since the phase 
boundary in Fig. 4 has a finite slope in the a-l plane, one 
can also induce the near-critical to supercritical transi- 
tion by increasing the /-value for a fixed a. 

For comparison, we also show in Fig. 4 the correspond- 
ing pha se boundary of the short-range model reported in 
Ref. [Tg] - As can be seen from the figure, the phase dia- 
gram of both the long-range and the short-range models 
are qualitatively similar. The near-critical phases in the 
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long-range model are replaced by the subcritical phases 
in the short-range model, and the phase boundaries of 
the long-range model tend to shift to larger values of a 
and to smaller values of I. 



B. THE MEAN DISPLACEMENT, THE MEAN 
NUMBER OF FAILED-BLOCKS AND THE 
MEAN STRESS-DROP 

The size of an earthquake event is usually measured 
by its magnitude. Other possible measures of event size 
might be the mean displacement Au, the mean number 
of failed-blocks AT& (corresponding to the size of rupture 
zone), and the mean stress-drop Af. In Figs. 5(a) and 
(b), we show the magnitude dependence of the mean dis- 
placement and of the mean number of failed-blocks for 
various values of a. An interesting observation is that 
the data in the near-critical regimes are grouped into 
two distinct branches, each corresponding to the small-a 
and large-a near-critical regions of Fig. 4. 

As can be seen from Fig. 5(a), the data in the small- 
er near-critical regime (a < a c \ ~ 2) lacks events of 
larger magnitudes and are characterized by smaller dis- 
placement, while those in the large-a near-critical regime 
(a > a C 2 ~ 25) are characterized by much larger dis- 
placement. All the data of the mean displacement Am in 
the near-critical regimes collapse, at least approximately, 
onto these two curves, which are both linear in the mag- 
nitude with a common slope ~ 0.01. Note that this slope 
is very small, indicating that the mean stress-drop in the 
near-critical regime hardly depends on the event magni- 
tude. This slope is also an order of magnitude smaller 
than the corresponding slope observed in the 2D short- 
range BK model, which was estimated to be about 0.1 

By contrast, the data in the supercritical regime {a c \ < 
<x < ot C 2) exhibit a significantly different behavior. At 
smaller magnitudes fi IS 5, they exhibit a crossover be- 
havior depending on its a-value between the two univer- 
sal near-critical curves: For smaller a close to a c %, the 
data tend to lie closer to the small-a near-critical curve, 
while for larger a close to a C 2, the data tend to lie closer 
to the large-a near-critical curve. The data in the super- 
critical regime suffer from significant finite-size effects at 
larger magnitudes /i c 5. The system-size N dependence 
of the data in the near-critical regime a = 30 is shown 
in the inset of Fig. 5 (a). As can be seen from the inset of 
Fig. 5(a), the data at larger magnitudes tend to level off 
as the system-size N is increased. 

The existence of the two near-critical curves and the 
crossover behavior between them are also clearly visible 
in Figs. 5(b) for the magnitude dependence of the mean 
number of failed-blocks iVj,. The two near-critical curves 
are again both strikingly linear with a common slope ~ 
0.99. The system-size N dependence is shown in the inset 
of Fig. 5(b) for the case of a = 30. As the system size is 
increased, the data at larger magnitudes tend to lie closer 



to a straight line of a slope ~ 0.99. 
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FIG. 5: The magnitude dependence of the mean displace- 
ment (a), the mean number of failed-blocks (b), and the mean 
stress-drop (c) of each seismic event of the 2D long-range BK 
model. In the main panels of Figs. (a) and (b), the frictional- 
parameter a is varied with fixing the system-size 60 x 60, 
while in the insets the system-size iV is varied for the case of 
a = 30. In Fig.(c), the system-size dependence of the mean 
stress-drop is shown for the case of a = 30. The parameters 
I and a are fixed to I = 3 and a — 0.01. 

In Fig. 5(c), we show the magnitude dependence of the 
mean stress-drop for the case of a = 30, with varying 
the system size N. Note that, although in the nearest- 
neighbor BK model the mean stress-drop of an event is 
essentially identical with (proportional to) the mean dis- 
placement of an event [l8| , such a simple relation between 
the mean displacement and the mean stress-drop does not 
hold in the present long-range BK model. Although a sig- 
nificant finite-size effect is observed, there clearly exists 
a tendency that the magnitude dependence becomes less 
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and less for larger systems. In real seismicity, the mean 
stress-drop is known to hardly depend on the event mag- 
nitude [l[. This suggests that the long-range nature of 
the elastic interaction of the crust might play a role in 
realizing the near-independence of the stress-drop on the 
event magnitude. We note that a similar independence 
was also observed in the mean-field-type ID long-range 
BK model studied by Xia et al 23j], and also in the ID 
long-range BK model with a power-law interaction stud- 
ied in Appendix B. 



THE LOCAL RECURRENCE-TIME 
DISTRIBUTION 
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FIG. 6: The log-log plot of the local recurrence-time distri- 
bution P(T) of large events of fi > fi c = 5 of the 2D long- 
range BK model for the frictional-parameter a — 10 and 30, 
each corresponding to the "supercritical" and "near-critical" 
regimes. The parameters I and a are fixed to I = 3 and 
a = 0.01. The system size is 160 x 80. The recurrence time T 
is normalized by its mean T, which is T = 31.5 and 9.98 for 
a = 10 and 30, respectively. The insets represent the semi- 
logarithmic plots including the tail part of the distribution. 
The tail part shows an exponential behavior for both cases of 
a = 10 and 30. 

In earthquake prediction, one natural quantity to be 
investigated might be the distribution law of the recur- 
rence time of large earthquakes. Characteristic earth- 
quake recurrence would mean the existence of character- 
istic time scales in earthquake recurrence, while critical 
earthquake recurrence would mean the absence of such 
characteristic time scales. Here, we study the nature of 
earthquake recurrence of the 2D long-range BK model 
via the local recurrence-time distribution function. 

In Figs. 6, we show on a log-log plot the computed local 
recurrence-time distribution function P{T) for the cases 
of a = 10 (a) and a = 30 (b), with fixing I = 3 and a = 
0.01. Each case corresponds to the supercritical and the 
near-critical regimes, respectively. The local recurrence 
time T is recorded when the next event occurs with its 



epicenter lying within distance r = 5 from the epicenter 
of the previous event. In the insets, the same data are re- 
plotted on a semi-logarithmic scale. The recurrence time 
is normalized by its mean T, which is Tv = 31.5 and 9.98 
for a = 10 and 30, respectively. As can be seen from the 
figure, P(T) exhibits an exponential tail at longer times 
for both cases of a = 10 and 30, with and without a peak 
structure at short times. 
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FIG. 7: The log-log plots of the local recurrence-time dis- 
tribution function P(T) of the 2D long-range BK model, 
with varying the magnitude-threshold /i c for the cases of the 
frictional-parameter a = 10 (a) and a = 30 (b). The in- 
sets represent the semi-logarithmic plots including the tail 
part of the distribution. The mean recurrence time is Tv = 
0.0262,14.3 and 31.5 (respectively for /i c = — oo,0 and 5) 
for a = 10, and f = 0.0135, 0.37 and 9.98 (respectively for 
/i c = — oo, and 3) for a = 30. 

In Figs. 7, we show P(T) for various values of the 
magnitude-threshold fi c , including the case of // c = — oo 
corresponding to no threshold at all (all events), for the 
cases of a = 10 (a) and of a = 30 (b). Both in the 
cases of a — 10 and 30, P{T) robustly exhibits a down- 
bending behavior for any choice of fi c . In the super- 
critical case of a = 10, a prominent peak observed at 
shorter T for \i c = 5 tends to be suppressed as \x c is taken 
smaller. In the near-critical case of a = 30, no character- 
istic peak is observed for any choice of \i c . The appear- 



8 



ance of the characteristic peak in P(T) at T ~ 0.01T 
in the supercritical regime and for larger events is well 
correlated with the appearance of the characteristic-peak 
component in the magnitude distribution R(fi) of Fig. 1. 
The recurrence-time distribution in real scismicity usu- 
ally does not exhibit a characteristic peak (see, e.g., Fig. 5 
of Ref.[30]). Hence, the behavior of P(T) in the near- 
critical regime seems closer to that of real seismicity. 
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FIG. 8: The log- log plot of the global recurrence-time dis- 
tribution function P(T) of the 2D short-range BK model. 
No constraint is imposed on the magnitude threshold, i.e., 
fi c = — oo. The system size is N x = 25 and N z = 100. The 
cases of a = 1.5 and 3.5, with l x = V3,l z — l,a = 0.01 
and with free boundary conditions applied both in x- and 
2-directions, precisely correspond to the cases studied in 
Ref . [241] . In these cases, we get Tv = 0.48 for a = 3.5, and 
Tv = 1.70 for a = 1.5. In the case of a = 3.5, the data 
taken under periodic boundary conditions applied in both al- 
and 2-directions as well as the data taken for the isotropic 
elastic couplings l w = l z = (1 + v / 3)/2 are also given. In the 
former case, we get T = 0.51, while in the latter case we get 
T — 0.48. The main panel represents the log-log plot of P(T), 
while the inset represents the semi-logarithmic plot including 
the tail part of the distribution. In all cases, the computed 
P(T) exhibits a down-bending behavior, in sharp contrast to 
Ref. [241] . The reason of this deviation is discussed in the text. 



BK model as studied by Hasumi, imposing no constraint 
on the distance between successive events. The result 
is shown in Fig. 8. First, we closely follow Ref.[24j by 
applying free boundary conditions in both directions on 
the lattice of size N x = 25 and N z = 100, assuming the 
anisotropic elastic parameters l x = \/?>, l z = 1, setting 
the other parameter values to a = 0.01 and a — 1.5 or 
3.5, and imposing no magnitude-constraint /i c = — oo. 
Precisely under these calculational conditions, Hasumi 
observed a critical straight-line P(T) for the case of 
a = 3.5, and a subcritical up-bending P(T) accompany- 
ing a characteristic larger-T peak for the case of a = 1.5. 
In sharp contrast to this, we observed here a subcritical 
down-bending P{T) for either value of a: See fig. 8. 

In the case of a = 3.5, we also examined the possible 
effect of the boundary conditions and of the anisotropy 
of the elastic constants on P(T) by simulating the 
model under periodic boundary conditions and with the 
isotropic elastic constants l x — l z — (l x + l z )/2, to find 
that the applied boundary conditions and the anisotropy 
of the elastic constants hardly affect the form of P(T) as 
shown in Fig. 8. 

In fact, Hasumi included the rise-time of earthquakes 
in his definition of the recurrence time, and assumed an 
extremely large loading rate, v — 10 -2 (25|. We believe 
that his choice of unrealistically large loading rate, com- 
bined with his definition of the recurrence time, is the 
cause of the deviation between our results and those of 
Ref. [24[ • As is well known, in real seismicity the loading 
rate is extremely small, being of order v = 10~ 8 — 10~ 9 . 
Then, with such a realistic choice of the rvalue, the 
recurrence-time distribution P(T) of the 2D BK model, 
either local or global, should behave in the way as re- 
ported in the present paper and in Ref. [IS], not as re- 
ported in Ref. [24j . irrespective of whether one includes 
the rise-time part of earthquakes in the definition of the 
recurrence time or not. 



D. TIME CORRELATIONS OF EVENTS 
ASSOCIATED WITH THE MAINSHOCK 



While the present results of P(T) turn out to be qual- 
itatively similar to those of the 2D short-range model 
calculated by the present authors [l8[, they differ sig- 
nificantly from the recent result of the recurrence-time 
(interoccurrence-time) distribution reported by Hasumi 
for the 2D short-range BK model (24[. Hasumi reported 
that the recurrence-time distribution P(T), defined glob- 
ally, exhibits either a supercritical, subcritical or critical 
behavior depending on the a-value, which is well cor- 
related with the behavior of the magnitude distribution 
function R{n). Such behaviors of P{T), however, were 
never observed in our calculation of the 2D BK model 
either in the short-range nor in the long-range case. 
In order to trace the cause of this significant deviation 
from Ref. [24[, we further calculated the global recurrence- 
time distribution on exactly the same 2D short-range 



In real seismicity, large events often accompany fore- 
shocks and aftershocks. In Fig. 9, we show the time cor- 
relation function between large events (mainshock) and 
events of arbitrary sizes, dominated in number by small 
events, for various values of the frictional-parameter a, 
with fixing / = 3 and a = 0.01. In the figure, we plot 
the mean number of events of arbitrary sizes occurring 
within 5 blocks from the epicenter of the mainshock be- 
fore (t < 0) and after (t > 0) the mainshock, where the 
occurrence of the mainshock is taken to be the origin of 
the time t = 0. The average is taken over all large events 
of their magnitudes of fi > fi c = 5. The number of events 
are counted here with the time bin of Atv = 0.02. 

As can be seen from Fig. 9, a remarkable acceleration of 
seismic activity occurs before the mainshock (t < 0) for 
a = 10 corresponding to the supercritical regime, while, 
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FIG. 9: The time correlation function of the 2D long-range 
BK model between large events of fi c = 5 (mainshock) occur- 
ring at time t = and events of arbitrary sizes, dominated in 
number by small events, occurring at time t for the cases of 
a = 10 and 30. The parameters I and <r are fixed to I = 3 and 
a = 0.01. Events of arbitrary sizes occurring within 5 blocks 
from the epicenter of the mainshock are counted. The neg- 
ative time t < represents the time before the mainshock, 
while the positive time t > represents the time after the 
mainshock. The average is taken over all large events with its 
magnitude fj, > fi c — 5. The system size is 160 x 80. 



for a = 30 corresponding to the near-critical regime, the 
time correlation is almost absent except for the suppres- 
sion of seismicity immediately before the mainshock. The 
behavior of the time-correlation function of the 2D long- 
range model turns out to be similar to those of the cor- 
responding 2D short-range model [HI . 



E. SPATIAL CORRELATIONS OF EVENTS 
BEFORE THE MAINSHOCK 
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FIG. 10: Event frequency before the large event of fi > [i c = 5 
(mainshock) plotted versus r, the distance from the epicenter 
of the upcoming mainshock, for several time periods before 
the mainshock in the 2D long-range BK model. The frictional- 
parameter a is a = 10 (a) and a = 30 (b) with / = 3 and 
a — 0.01, each corresponding to the "supercritical" and "near- 
critical" regimes. The system size is 160 x 80. The insets 
represent similar plots at shorter times. 



In this subsection, we examine the time- development 
of spatial seismic correlations before the mainshock of 
fi > [i c = 5. In Figs. 10, we show the spatial seismic cor- 
relation functions between the mainshock and the pre- 
ceding events of arbitrary size, dominated in number by 
small events, for several time periods before the main- 
shock, with fixing I — 3 and a = 0.01. Figs. 10(a) and (b) 
represent the cases of a = 10 in the supercritical regime 
and of a = 30 in the near-critical regime, respectively. 
Insets represent shorter-time behaviors. 

As can be seen from Fig. 10(a), for a = 10, the fre- 
quency of small events are enhanced preceding the main- 
shock at and around the epicenter of the upcoming main- 
shock. For small enough t, such a cluster of smaller 
events correlated with the large event may be regarded as 
foreshocks. Just before the mainshock, the frequency of 
smaller events is suppressed in a close vicinity of the up- 
coming mainshock, while it continues to be enhanced in 
the surrounding blocks, a phenomenon closely resembling 



the "Mogi doughnut" HSHll]. The spatial range where 
the quiescence occurs is narrow, only of a few blocks. 

For a = 30, as can be seen from Fig. 10(b), the seis- 
mic acceleration preceding the mainshock is hardly dis- 
cernible, while the doughnut-like quiescence is still real- 
ized. We note that the quiescence just before the main- 
shock is robustly observed in the BK model, independent 
of its dimensionality, the interaction range and the pa- 
rameter values. Indeed, it has been observed both in ID 
and 2D, both with the short-range and long-range inter- 
actions [reLlrTj. 



F. SPATIAL CORRELATIONS OF EVENTS 
AFTER THE MAINSHOCK 

In this subsection, we examine the time-development 
of spatial seismic correlations after the mainshock of /j, > 
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fj, c — 5. Figs. 11(a) and (b) represent the cases of a = 10 
in the supercritical regime and of a = 30 in the near- 
critical regime. Insets represent shorter-time behaviors. 
Computational conditions are taken to be the same as in 
Figs.10. 
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FIG. 12: The local magnitude distribution of the 2D long- 
range BK model for several time periods before the mainshock 
of /i > /j, c = 3 for the cases of a — 10 (a) and a = 30 (b), 
each corresponding to the "supercritical" and "near-critical" 
regimes respectively. Events whose epicenter lies within 5 
blocks from the epicenter of the upcoming mainshock are 
counted. The parameters I and a are fixed to I = 3 and 
a = 0.01. The system size is 160 x 80. In (a), an apparent 
B-value decreases before the mainshock, while in (b) it stays 
almost unchanged. 



FIG. 11: Event frequency after the large event of /i > [i c — 5 
(mainshock) plotted versus r, the distance from the epicenter 
of the preceding mainshock, for several time periods after the 
mainshock in the 2D long-range BK model. The frictional- 
parameter a is a = 10 (a), and a = 30 (b) with I = 3 and 
a — 0.01, each corresponding to the "supercritical" and "near- 
critical" regimes. The system size is 160 x 80. The insets 
represent similar plots at shorter times. 



As can be seen from the figures, spatiotemporal seis- 
mic correlations are almost absent after the mainshock 
in both cases of a = 10 and a — 30. In the present 
2D long-range model, the event frequency hardly changes 
with distance r nor with time t even in the supercritical 
regime, in contrast to the short-range case where non- 
trivial spatiotemporal correlations are observed to some 
extent even after the mainshock Hal. 



G. THE TIME-DEPENDENT MAGNITUDE 
DISTRIBUTION BEFORE THE MAINSHOCK 

In real seismicity, an appreciable change of the B- value 
of the magnitude distribution has been reported preced- 
ing larg e earthquakes: Often a decrease of the B-value 
[H, [29l[30j , but sometimes an increase of it [3l[ . Obvi- 
ously, a possible change in the magnitude distribution 
preceding the mainshock possesses a potential impor- 
tance in earthquake prediction. 

In Figs. 12, we show the "time- resolved" local mag- 
nitude distributions for several time periods before the 
mainshock for the cases of a = 10 (a) and of a = 30 
(b), with fixing / = 3 and a = 0.01. Only events with 
their epicenters lying within 5 blocks from the upcoming 
mainshock of /i > fi c = 3 are counted here. 
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As can be seen from Fig. 12(a), in the supercritical 
regime, an apparent B- value describing the smaller mag- 
nitude region, fj, < 2, gets smaller as the mainshock is 
approached. Indeed, the B-value is reduced from the all- 
time value B ~ 1.33 to B ~ 1.11 here. By contrast, as 
can be seen from Fig. 12(b), an apparent B-value hardly 
changes in the near-critical regime even when the main- 
shock is approached. It stays at around B ~ 0.59. 

Likewise, one can also study the "time-resolved" local 
magnitude distributions after the large event. Seismic 
events are quite scarce after the large event, however, 
as can be seen from Figs. 9 and 11, i.e., few aftershocks 
observed. Hence, it is statistically difficult to obtain the 
"time-resolved" local magnitude distributions after the 
mainshock in the present model. 



IV. SUMMARY AND DISCUSSION 

Spatiotemporal correlations of the 2D BK model with 
the long-range interaction were studied by means of ex- 
tensive numerical computer simulations. The long-range 
interaction, which takes account of the effect of the elastic 
body adjacent to the fault plane, falls off with distance r 
as 1/r 3 . Seismic properties of the model can be summa- 
rized in the form of a phase diagram shown in Fig. 4. The 
2D long-range model turns out to possess three distinct 
"phases", i.e., a small-a near-critical phase, a supercrit- 
ical phase and a large-a near-critical phase. 

The long-range BK model is certainly a more faith- 
ful representation of an earthquake fault than the short- 
range BK model. Although some of the properties are 
more or less common between in the short-range and in 
the long-range models, several important differences exist 
in several observables. 

Generally speaking, spatial seismic correlations of the 
long-range models tend to be suppressed compared with 
those of the short-range models. Such a suppression of 
spatial seismic correlations in the long-range model might 
intuitively be understandable, because the long-range na- 
ture of the interaction serves to smear out the spatial 
variation. 

Most interestingly, it has been found that the magni- 
tude distribution of the 2D long-range model exhibits 
a "near-critical" power-law-like behavior close to the 
Gutcnberg-Richter law, for a wide parameter range with 
its -B-value insensitive to the model parameter, B ~ 0.55, 
in sharp contrast to the cases of the 2D short-range 
model and of the ID short-range and long-range mod- 
els where such a near-critical behavior is realized only by 
fine-tuning the model parameter to a special value. Since 
the GR law is known to be robustly observed over differ- 
ent fault zones of varying locations and depths, possibly 
characterized by varying material parameters, the power- 
law feature of the magnitude distribution should be a 
stable attribute of earthquake occurrence, not a special 
property requiring a fine-tuning of the material parame- 
ter. In that sense, stable occurrence of the near-critical 



magnitude distribution over a wide parameter range in 
the 2D long-range BK model might be of relevance to real 
seismicity. The observed B-value, B ~ 0.55 (b ~ 0.83), 
is not far from the one observed in real faults B ~ 2/3 
(b ~ 1). It should be noticed that the near-critical magni- 
tude distribution observed here is not a truly critical one, 
since, for sufficiently large magnitude, the magnitude dis- 
tribution falls off sharply. The apparent power-law-likc 
behavior does not extend toward larger magnitudes in- 
definitely. 

In real seismicity, there holds an empirical law that 
the mean stress-drop of an earthquake is nearly constant 
irrespective of the event magnitude. Although in the 
short-range BK models, the mean stress-drop increases 
considerably as the magnitude gets larger, in the long- 
range BK models the mean stress-drop hardly depends on 
the event magnitude in a wide parameter range. Hence, 
the long-range BK model in the near-critical regime has 
an obvious advantage that it can reproduce the observed 
constancy of the stress-drop. 

Large events of the long-range model usually accom- 
pany foreshocks together with a doughnut-like quiescence 
as their precursors, while they hardly accompany after- 
shocks with almost negligible seismic correlations ob- 
served after the mainshock. Such absence of post-seismic 
activity correlated with the mainshock is more prominent 
in the long-range model than in the short-range model. 
Concerning pre-seismic activity preceding the mainshock, 
an appreciable change of the effective _B-value has occa- 
sionally been observed both in the long-range and short- 
range models in 2D. The B-valuc is cither increased, de- 
creased or unchanged in 2D, depending on the system is 
in the subcritical, supercritical or near-critical phase. 

In this way, the long-range 2D BK model, which is ap- 
parently the most realistic version among the types of 
the BK model studied so far, appear to give a reasonable 
description of real seismicity in its near-critical regime. 
The model can explain the GR-like magnitude distribu- 
tion with the B-value, B ~ 0.55, stably realized over a 
rather wide parameter range, the near independency of 
the stress drop on the event magnitude and the absence 
of a characteristic peak in the recurrence-time distribu- 
tion of earthquakes, etc. Meanwhile, characteristic fea- 
tures become most eminent in the supercritical regime, 
particularly for large events. 

Not all properties of real seismicity, however, are ex- 
plained by the model. For example, the Omori law fre- 
quently observed in real seismicity cannot be reproduced 
in the BK model. This may suggest that the effects not 
taken account in the present model, e.g., processes like 
the water migration through the crack, the slow chemi- 
cal process at the fault or the elastoplasticity associated 
with the ascenosphere, are important in realizing the af- 
tershock obeying the Omori law. 

In order to make a further link between the BK model 
and the real world, we estimate here various time and 
length scales involved in the BK model. For this, we 
need to estimate the units of time and length of the BK 
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model in terms of real-world earthquakes. Concerning 
the time unit w -1 , we estimate it via the rise time of large 
earthquakes, ~ tt/ui, which is typically about 10 seconds. 
This gives an estimate of w _1 ~ 3 sec. Concerning the 
length unit C, we estimate it making use of the fact that 
the typical displacement in large events of our simulation 
is of order one C unit, which in real- world large earth- 
quakes is typically 5 meters. Then, we get C ~ 5 meters. 
Since the loading rate v' associated with the real plate 
motion is typically 5 cm/year, the dimensionless loading 
rate v = v'/(£ui) is estimated to be v ~ 10~ 9 . 

In our simulation of the BK model, the doughnut-like 
quiescence was observed before the mainshock at the time 
scale of, say, tv IS 10~ 2 — 10 _1 . This time scale cor- 
responds to about 1-10 years. In our simulation, the 
doughnut-like quiescence was observed in the region only 
within a few blocks from the epicenter of the mainshock. 
To give the corresponding real- world estimate, we need 
the real- world estimate of our block size a'. In the BK 
model, the length scale a' is entirely independent of the 
length scale £, and has to be determined independently. 
We estimate a' via the typical velocity of the rupture 
propagation, la'cu, which is about 3 km/sec in real earth- 
quakes. From this relation, we get a 1 ~ 3 km. The length 
scale associated with the doughnut-like quiescence is then 
estimated to be 3 ~ 6 km in radius. 

Throughout our present simulations, we have assumed 
the velocity-weakening friction force. The extent of the 
velocity weakening is mainly described by the parameter 
a, which has a dimension of the inverse velocity. The 
unit of a^ 1 is then estimated to be ~ 1 m/sec. The 
recent high-velocity friction measurements indicate that 
the friction coefficient of serpentinite drops significantly 
at the slip velocity of order 0.1 ~ 1 m/sec [321 ] ■ This 
roughly corresponds to the a- value of order 1 a S; 10, 
which is indeed the values of interest here. 

The seismic moment M is approximately given by, 

M ~ GDS, (9) 

where G is the shear modulus, D is the displacement and 
S is the rupture-zone size. The shear modulus G of the 
crust is typically about 30 GPa. The moment magnitude 
to measured in the MKS unit is defined by, 

to = H l ogl0 M - 6. (10) 

From Eq.(7.1) and (7.2), we can obtain the relation be- 



tween the moment magnitude to in real world and the 
magnitude // in the BK model as, 

to ~ 0.29^ + 6.09. (11) 

In the BK model, an upper cut-off magnitude [i m ax of 
the GR-like behavior is found to be about 6 5 Umax ~ 8, 
which corresponds to 7.8 IS m max 5 8.4 in real world. Al- 
though it is tempting to speculate that the "interrupted 
power-law" or "near-criticality" observed in our model 
simulation might somehow be related to real observa- 
tion, it is not known whether there really exists an upper 
cut-off magnitude in real seismicity. 

In the long-range BK model, the mean stress-drop 
hardly depends on the event magnitude in the near- 
critical regime. In the 2D long-range BK model in the 
near-critical regime, the mean stress-drop was estimated 
to be At ~ 1.6. This corresponds in the real world the 
mean stress drop of ~ 80 MPa. Since the mean stress- 
drop is 1 ~ 10 MPa in real seismicity, the estimated value 
of the mean stress-drop is a bit larger than but roughly 
consistent with the real value. 

The present study was performed under many assump- 
tions, e.g., an earthquake fault is completely flat, ma- 
terial parameters are homogeneous, there is no depth 
dependence in the material parameters, a friction force 
depends on the velocity alone, etc. As one of such as- 
sumptions, we have employed a static approximation in 
our simulation of the long-range BK model, i.e., we have 
assumed that the velocity of the seismic- wave propaga- 
tion is sufficiently larger than the rupture velocity. In real 
earthquakes, however, the rupture velocity is comparable 
to the shear-wave velocity. Thus, it is clearly desirable 
to perform simulations based on a fully dynamical elastic 
theory 

In spite of such limitations of the model, our present 
study has revealed that the 2D long-range BK model can 
reproduce several important aspects of real seismicity. 
We hope that the present analysis might give a step to- 
ward the fuller understanding of the statistical properties 
of earthquakes. 
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APPENDIX A: THE DERIVATION OF THE 
INTERACTION BETWEEN TWO ARBITRARY 
BLOCKS 



where U n (x) represents a displacement in the n-th di- 
rection at a spatial point x=(x\, x 2 , x 3 ) in the elastic 
body, AUi(x') is a relative displacement in the i-th di- 
rection across the fault surface E, rij(x') is the normal 
unit vector on the fault surface, and C'ij pq is an clastic 
constant. The Green's function G np (x;x') is a displace- 
ment in the n-th direction at a point x—(x±, x 2 , X3) due 
to a unit force acting along the p-th direction at a spatial 
point x' = (x[, x' 2: x 3 ). We assume the fault surface to be 
the xia^-plane which slips only in the xi-direction. The 
elastic body is assumed to be isotropic, homogeneous and 
infinite. 

We consider a static version of the Navier's equation 
as a differential equation describing the elastic body, 



(A+ u )— ( dGjn 
dxi \ dx* 



d fdGi, 



dxj \ dxj 



-S in 5(xi)S(x 2 )S(x 3 ), 
(A2) 



the associated Green's function being given by 

A + fi d 2 R 



6 n p 1 



Airfi R 8"7r^(A + 2fi) dx n dx p 



(A3) 



R=\x — x'\ = yj(xi ~ <) 2 + (x 2 - 4) 2 + (a* - 4) 2 , 

where d np is the Kronecker's delta, and A and fi are 
Lame's constants. 

The stress tensor nj is related to the strain tensor 
via the Hooke's law, 



Aefefei5y + 2/xe, 



2 \ dxi dxk 



(A5) 



(A6) 



Then, we consider the situation where an infinitesimal 
part of the fault plane dT,' located at (x' 1 ,0,x' 3 ) slips by 
an amount AU. By using Eqs.(Al), (A3), one gets the 
stress on the x 1X3 -plane (x 2 = x' 2 ) as 



n 2 (xi - x[ 1 0,x 3 - x' 3 ) = 



/j(3A+2^)dE 
47r(A+2^) 



Ro — 



(si-Si) 1 1ji (S3-S3) 



AC/, 



yj \x 1 - x[y + {x 3 - x' 3 y 



(A7) 



(A8) 



In this appendix, based on an elastic theory, we derive 
the effective interaction between two arbitrary blocks on 
a fault plane in the 2D BK model. We begin with the 
static version of the representation theorem [331 ] that rep- 
resents a displacement in an elastic body induced by a 
slip on a crack surface (or a fault plane) as 

U n (x)= J J ^Ui{x')C ijm nj(x')-^-G np {x-,x') dE', 

(Al) 



The result means that the stress decays with distance Ro 
as 1/Rq on the fault plane. Indeed, Maruyama discussed 
a static version of the three-dimensional source mechanics 
of earthquakes 34]. The result we have obtained here 
corresponds to his result. 

Now, we wish to apply Eqs. (A7) and (A8) derived 
from an elastic theory to the BK model. First, we dis- 
cretize the fault plane into blocks of linear size a'. Sec- 
ond, we regard AU to be a relative displacement between 



two blocks, i.e., AU = Ui. 



where Uij denotes a 
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displacement of a block at a site (i, j). Note that, by this 
choice of At/, one has a vanishing self-interaction, since 
the relative displacement with itself always vanishes. 

The spring constant K(i, , j')(Uij — Ui>j>) = a' 2 r 
between the blocks at a site and at (i',f) is then 

given by 

K(i,j;i',j') = C^-P^ + l' z 2( 2-p^, (A9) 



7/2 



fj,(3X + 2fj,)a' 
4tt(A + 2^) ' 



2/xV 
4tt(A + 2^)' 



(A10) 
(All) 



Then, after the block discretization and the replace- 
ment At/ = Ui — Ui>, the spring constant defined by 
K(i;i')(Ui — Ui>) = cl't is obtained as 



K(i;i') = r 



1 



l\2 ■ 



/2 _ fj,(\ + fj,)a' 
~ tt(A + 2/z) ' 



(A20) 



(A21) 



The dimensionlcss inter-block interaction is then given 



by 



\i — I'r 



(A22) 



(A12) 



by 



The dimensionlcss inter-block interaction is then given 



2 {i-i>? 2 (i-jQ 2 



){ui,j - Ui>j>), (A13) 



/ 2 — JL. — 

<>x — , — 



li(3X + 2n)a' 



k p 47r(A + 2^)fc P 



I 2 - 



2/iV 



47r(A + 2^)fc p 



(A14) 



(A15) 



where k p is the spring constant introduced in §3, and Uij 
is a dimensionless displacement defined in §3. 

Similarly, for the long-range ID BK model, we can ob- 
tain the interaction K(i; i') between two arbitrary blocks 
at site i and i' . In the long-range ID BK model, we have 
assumed the fault and the elastic body to be a rigid body 
in the ^-direction. In this case, a static version of the 
Navier's equation may be written as 



(A + /x) 



d fdGi 



d fdd 



-S in d{x 1 )S(x 2 ), 
(A16) 



' dxi \ dxj J ' dxj \ dxj 

the associated Green's function being given by 

A + ^ <9 2 [i? 2 (logi?- 1)] 



Gn P (x;x') = -^lo g R+- 8MX + 2 ^ 



(A17) 



R 



\x-x'\ = ^(x 1 -x' 1 ) 2 + (x 2 -x> 2 )2. (A18) 



By using Eqs.(Al) and (A13), one obtains the stress on 
the xia;3-planc (x 2 — x' 2 ) as 



t 12 (x 1 - x[,0) 



H(X + fj,)dE 



1 



ir{X + 2fi) (xi~x[)' 



: AU. (A19) 



r _ n(X + n)a' 
kp Tr(X + 2fi)k p ' 



(A23) 



APPENDIX B: THE ID BK MODEL WITH THE 
LONG-RANGE INTERACTION 

In this appendix, we show some of the results of our 
numerical simulations on the ID BK model with the long- 
range interaction decaying as 1/r 2 . 

Typical behaviors of the magnitude distribution are 
shown in Figs.l for the case of I = 3 and a — 0.01. 
Figs. 13 (a) and (b) exhibit for smaller and larger 

a, respectively. The peculiarity of the ID long-range 
BK model is that, for sufficiently small values of a 5 
0.6, only one-block events occur under periodic boundary 
condition in the steady state realized after transients. 
Under free boundary condition, on the other hand, such 
an exclusive occurrence of one-block events does not arise 
for any a. We note here that the behavior for smaller a 
is rather sensitive to the choice of the time discretization 
At. In the region of smaller a, we need to take At as 
small as 10~ 5 to get stable results. Otherwise, totally 
different behaviors would sometimes arise. 

In the range of 0.7 IS a IS 1, events involving more 
than one block begin to occur, where the associated R(fj) 
exhibits a "subcritical" behavior bending down rapidly 
at larger magnitudes, as can be seen from Fig. 13(a). As 
a is increased, weights of larger events tend to increase 
gradually, and at a = 1, R((a) exhibits a near straight-line 
"near-critical" behavior close to the GR-law behavior. 

As a is increased further beyond a = 1, i?(yu) devel- 
ops a characteristic peak and exhibits a "supercritical" 
behavior, deviating from the GR law at larger magni- 
tudes fj, ^ fb ~ 1, while it still exhibits a near straight- 
line behavior corresponding to the GR law at smaller 
magnitudes /j, IS p.- As a is further increased, the peak 
at a larger magnitude becomes less pronounced, and at 
a ~ 15, R{p.) exhibits a near-critical behavior again with- 
out a characteristic peak. For a c 15, R(/J.) exhibits a 
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FIG. 13: The magnitude distribution of earthquake 

events of the ID long-range BK model for the parameters 
I — 3 and a — 0.01. Fig. (a) represents R(fJ,) for smaller values 
of the frictional parameter < a < 5, while Fig.(b) represents 
R(n) for larger values of the frictional parameter 5 < a < oo. 
The system size is N = 800. 



subcritical behavior, rapidly bending down at larger mag- 
nitudes. Finally, events involving more than one block 
suddenly disappear. In the range of a ^ 24, only one- 
block events occur. As in the case of smaller a, we need 
to take the time discretization At sufficiently small in or- 
der to correctly reproduce such a behavior in this regime 
of larger a. 

While the magnitude distributions presented here are 
the first data on the ID BK model with the 1/r 3 long- 
range interaction, we wish to make some comparison with 
the earlier data for the related ID BK models. The 
magnitude distribution of the ID short-range (nearest- 
neighbor) BK model was studied by several authors, in- 
cluding the earlier calculation of Carlson, Langer and 
collaborators Isjjly] as well as of our own [18]. The 
data of Ref.0, Q corresponded to the "supercritical" 
regime (a — 2.5,3 and 4) and the "near-critical" regime 
(a = 1). Our present data are qualitative similar to 
those of Refs.fEl, [g] in these regimes, though the GR-like 
behavior at smaller magnitudes, i.e., the linearity of the 



R{p) curve, seems less pronounced in our present case 
and in Ref.pJ] than in Ref.pl [f|. This is due to the dif- 
ferent choice of the Z-value: Carlson et al took I to be 
large 6 ~ 14, while we mostly choose I — 3 here and in 
Ref.[Hj]. 

By contrast, if we compare our present for 
the 1/r 3 long-range BK model with the one obtained 
in Ref.[22j, |23j for the mean- field- type long-range BK 
model, there exists some appreciable qualitative differ- 
ence. Namely, even in the "supercritical" regime of a = 2 
and 2.5, the magnitude distribution of Ref . [22l. |23| ex- 
hibits no characteristic peak at a larger magnitude, but 
rather exhibits a down-bending "subcritical" -type behav- 
ior. In Ref . [H, H3| , a characteristic peak in R(/j.) is dis- 
cernible in the region of smaller a (a = 0, 0.5) where 
we have observed here either "one-block events only" be- 
havior or "subcritical" behavior without a characteristic 
peak. We have checked that this qualitative difference is 
not due to the different choice of the /-value in the two 
calculations. Thus, the behavior of the magnitude dis- 
tribution appears to differ substantially between in the 
mean- field- type long-range model and in the 1/r 3 long- 
range model. 




a 

FIG. 14: The phase diagram of the ID BK models with 
the long-range interaction in the frictional-parameter a versus 
elastic-parameter I plane. The parameter a is set to a = 0.01. 
To draw a phase diagram, the parameter range < a < oo 
and 1 < Z < 10 is studied by simulations. 

In Fig. 14, we summarize the behavior of R(n) of the ID 
long-range BK model in the form of a "phase diagram" in 
the frictional-parameter a versus the elastic-parameter / 
plane for the case of a = 0.01. The phase diagram con- 
sists of five distinct regimes, two of which are "one-block 
events" regimes, two are "subcritical" regimes and one 
is a "supercritical" regime. The transition between the 
small-a subcritical regime and the supercritical regime 
appears to be continuous (gradual), in contrast to the 
one of the 2D long-range model. The transition between 
different "phases" is primarily dictated by the a-value. 
Since the "phase boundary" in Fig. 14 has a finite slope 
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in the a-l plane, one can also induce the transition by 
increasing the lvalue for a fixed a. 

In the mains panels of Figs.l5(a)-(c), we show the mag- 
nitude dependence of the mean displacement, the mean 
number of failcd-blocks and the mean stress-drop of the 
ID long-range BK model for various values of a. In the 
insets, we show the system-size dependence of each quan- 
tity for the case of a = 1. 
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FIG. 15: The magnitude dependence of the mean displace- 
ment (a), the mean number of failed-blocks (b), and the mean 
stress-drop (c) of each seismic event of the ID long-range BK 
model. In the main panels, the frictional-parameter a is var- 
ied with fixing the system-size N = 800, while in the insets 
the system-size N is varied for a = 1. The parameters I and 
a are fixed to I = 3 and a = 0.01. 

As can be seen from Figs. 15(a) and (b), the data might 



roughly be grouped into three different categories, each 
corresponding to the small-a subcritical regime, the su- 
percritical regime and the large-a subcritical regime, al- 
though the transition between these behaviors is rather 
gradual. As compared with the corresponding 2D mod- 
els, including both the short-range model studied in [Lsl ] 
and the long-range model studied in §3, the scaling prop- 
erty is much more obscured here in ID. The data in the 
subcritical regimes do not collapse on top of each other, 
nor exhibit a straight-line power-law-like behavior. 

As can be seen from Fig. 15(c), the mean stress-drop of 
a seismic event Af hardly depends on its magnitude /i 
except for large earthquakes. There is even a tendency 
that the mean stress-drop becomes more independent of 
the event magnitude as one studies larger systems (see 
the inset). Similar independence is also observed in the 
2D long-range model in §3, as well as in Ref.[23| for the 
mean-field-type ID long-range model, and might be con- 
trasted to the property of the corresponding short-range 
model where the mean stress-drop exhibits more pro- 
nounced magnitude dependence [181 ] . 
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FIG. 16: The local recurrence-time distribution function 
P(T) of the ID long-range BK model of the frictional pa- 
rameter a — 1, with varying the magnitude threshold /i c . 
The main panels represent the log-log plot of P(T) and the 
insets represent the semi-logarithmic plots including the tail 
part of the distribution. The mean recurrence time T is 
T = 0.000916,0.369 and 19.7, respectively for ju c = -oo,0 
and 3. 

We have also computed the local recurrence-time dis- 
tribution P(T) for events of their magnitude /i > // c . The 
local recurrence time T is defined by the time passed until 
the next event occurs with its epicenter lying in a vicin- 
ity of the previous event within distance of r = 30-blocks 
from the epicenter of the previous event. The behavior of 
the computed local recurrence-time distribution is qual- 
itatively similar to the one of the 2D long-range model 
given in §3; an exponential tail at longer T, with or with- 
out a characteristic peak at shorter T in the supercritical 
or in the subcritical regimes, respectively. 

One case of interest in the ID model might be the 
a = 1 near-critical case located at the phase boundary 
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between the small-a subcritical regime and the super- 
critical regime, since such a region is absent in the cor- 
responding 2D model due to the discontinuous nature of 
the transition. Thus, in Fig. 16, we show on a log-log plot 
the computed P(T) for the case of a = 1, with fixing 
I = 3 and a = 0.01, for various values of the magnitude 
threshold fi c . As can be seen from the figure, P(T) tends 
to exhibit a power-law-like behavior at larger T as the 
magnitude threshold /x c is taken smaller. The associated 
exponent is estimated to be about ~ 2.6. This suggests 
that, at a = 1, the occurrence of small events has a crit- 
ical feature, while such a critical feature is weakened for 
larger events. Such a critical feature was not seen in the 
recurrence-time distribution of the 2D long-range model 
studied in §3. We note that, even in ID, such a critical 
P(T) is realized only at a = 1. For other values of a, 
P(T) robustly exhibits an exponential tail at longer T 
(not shown here). 

We have also calculated various spatiotemporal corre- 
lation functions for the ID long-range BK model, most of 
which show behaviors qualitatively similar to the ones ob- 
served for the 2D long-range BK model. Among them, we 
show in Figs. 17 the "time-resolved" local magnitude dis- 
tributions for several time periods before the large event 



for the cases of a = 1 (a), a = 5 (b) and a = 15 (c), with 
fixing / = 3 and a = 0.01. Only events with their epicen- 
ters lying within 30 blocks from the upcoming mainshock 
is counted here. We define the mainshock as a large event 
of /i > /i c = 3. 

For the case of a = 1, as shown in Figl7(a), an ap- 
parent B-value describing the smaller magnitude region 
fj, S — 1 gets smaller from the all-time value B ~ 0.79 
to the short-time value B ~ 0.67 as the mainshock is 
approached. Such a decrease of the -B-value is opposite 
to the one observed in the corresponding ID short-range 
model at a = 1 where the B-value gets larger as the 
mainshock is approached [H, [llj • 

For the case of a = 15, by contrast, an apparent B- 
value describing the smaller magnitude region gets larger 
as the mainshock is approached: See Fig. 17(c). 

For the case of a = 5, the time development of the 
magnitude distribution exhibits a somewhat different be- 
havior as shown in Fig. 17(b). As the mainshock is 
approached, the magnitude distribution R(fi) is devel- 
oped from the supercritical all-time behavior to the near- 
critical straight-line behavior characterized by a slope 
5 ~ 1.11. 
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FIG. 17: The local magnitude distribution of the ID long- 
range BK model for several time periods before the mainshock 
of jit > fi c = 3, for the cases of a = 1 (a), a = 5 (b), and 
a = 15 (c). Events whose epicenter lies within 30 blocks 
from the epicenter of the upcoming mainshock are counted. 
The parameters I and a are fixed to I — 3 and a — 0.01. The 
system size is iV = 800. In (a), the apparent B- value decreases 
before the mainshock, while, in (c), it increases before the 
mainshock. 



